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Quantum Monte Carlo methods are used to study a quan- 
tum phase transition in a ID Hubbard model with a stag- 
gered ionic potential (A). Using recently formulated methods, 
the electronic polarization and localization are determined 
directly from the correlated ground state wavefunction and 
compared to results of previous work using exact diagonal- 
ization and Hartree-Fock. We find that the model undergoes 
a thermodynamic transition from a band insulator (BI) to a 
broken-symmetry bond ordered (BO) phase as the ratio of 
(7/ A is increased. Since it is known that at A = the usual 
Hubbard model is a Mott insulator (MI) with no long-range 
order, we have searched for a second transition to this state 
by (i) increasing U at fixed A and (ii) decreasing A at fixed 
U. We find no transition from the BO to MI state, and we 
propose that the MI state in ID is unstable to bond ordering 
under the addition of any finite ionic potential A. In real 
ID systems the symmetric MI phase is never stable and the 
transition is from a symmetric BI phase to a dimerized BO 
phase, with a metallic point at the transition. 

I. INTRODUCTION 

Strongly-correlated systems of interacting electrons 
lead to many of the most interesting phenomena observed 
in solid state physicsEI. As a function of the interaction 
strength, there can be quantum phase transitionsS char- 
acterized by an order parameter with the possible devel- 
opment of long-range order and a transition to a broken 
symmetry state. Interactions can also lead to "Mott in- 
sulators" (MI) and to metal-insulator transitionsi. An 
important question is whether or not in the thermody- 
namic limit a Mott insulator must be associated with a 
phase transition that is accompanied by a broken sym- 



metry and a corresponding order parameter. In his orig- 
inal work, Motti argued that the insulating character 
did not depend upon an order parameter. On the other 
hand, SlaterS emphasized the relation of the insulating 
behavior to the long range order, and in many cases it 
is known that the MI state must be accompanied by a 
broken symmetry!. 

To address such issues theoretically we must have 
methods that can clearly distinguish metals from insula- 
tors, i.e., the ability to transport chargei~i vs. localiza- 
tion of the electronJi. Insulators at absolute zero can not 
transport arbitrary amounts of charge macroscopic dis- 
tances across their bulk; however, the center of electronic 
charge can shift in response to external fields, which is 
described in terms of changes in polarizatiorBS. The po- 
larizability is characterized by the degree of electronic 
delocalizationil which increases with the proximity to the 
metallic state. Recently, there have been new develop- 
ments defining macroscopic polarization and localization 
in terms of the insulating ground state wavefunction^^lli. 
These theories formulate the polarization and localiza- 
tion in terms of Berry's phasesE^I which can be calcu- 
lated using "twisted boundary conditions" or in terms 
of the expectation value of an exponentiated operator. 
Such twisted boundary conditions have been applied in 
the past to study metals and approach metal-insulator 
transitions from the metallic side0~@'i. With the re- 
cently developed methods for insulators, there are now 
complementary toolfS to provide quantitative informa- 
tion on the divergence of the localization length as one 
approaches the metal-insulator transition from the insu- 
lating side. 

Generalized Hubbard models00 are well-suited for 
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studies of fundamental issues regarding metals and insu- 
lators because they are simple models that exhibit a wide 
range of behaviors depending upon the parameters in the 
models. The simplest of all, the original Hubbard model 
with an on-site interaction U and nearest neighbor hop- 
ping < in ID, was solved exactly by Licb and WuEl. Their 
paper, entitled "Absence of Mott Transition conveys 
the point that there is no change of spatial symmetry and 
no phase transition at any positive U. At half-filling the 
model is metallic at U — 0, whereas at any positive in- 
teraction U a gap exists to charge excitations but no gap 
exists to spin excitations. This is commonly referred to 
as the MI state, but in this case there is no "Mott transi- 
tion" . At any other filling, the model is always metallic. 
There is never a state that would be called an ordinary 
band insulating (BI) state. However, in systems of higher 
dimensionality (d > 2) , a MI state is always accompanied 
by a broken symmetry!^. 

Many new possibilities emerge for generalized Hubbard 
models in ID. The ionic ID Hubbard model with two 
inequivalent sites, proposed by NagaosaEl and later by 
Egami0 as a model ferroelectric, is ideal for studying 
how quantized particle transport is modified by electron 
correlation in a many body system. On general grounds 
we expect a transition to occur from an ionic band in- 
sulator to a strongly correlated Mott insulator as U is 
increased. Evidence for such a transition was found in 
exact-diagonalization calculationsEi'0, where the elec- 
tronic polarization was found to jump abruptly between 
two discrete values fixed by the existence of two centers 
of inversion at the two sites. Such behavior has been 
termed a "topological transition''^ that occurs in finite 
systems and therefore is distinct from a true quantum 
phase transition. These solutions predict that the model 
has a metallic point separating two insulating phases and 
that a ferroelectric polarization results only if the atomic 
sites are displaced from the centers of inversion. 

However, recently Fabrizio, et al.,!! have proposed that 
this model will instead exhibit two quantum phase tran- 



sitions: one from a BI state to a long range bond ordered 
(BO) state, predicted to be in the Ising universality class, 
and a second from the BO to the MI state, predicted to be 
a Kosterlitz-Thouless transition. Such transitions to BO 
states have recently been found in ID Hubbard models 
with extended interactions (U-V) by Nakamura.Ei The 
BO state is a broken symmetry state in which the sys- 
tem becomes ferroelectric due strictly to electron-electron 
interactions even if all the atoms are at centers of inver- 
sion. 

During the course of the present work, two preprints 
have reported calculations of charge and spin gaps in the 
mode]ii@. Even though each work uses the density ma- 
trix renormalization group (DRMG) that allows studies 
of very large ID systems, each group reports great dif- 
ficulty in extrapolating to large size the small spin gaps 
and the two papers come to opposite conclusions regard- 
ing the existence of the BO state. 

The purpose of this paper is to study the ionic Hubbard 
model using a method that (i) will treat electron correla- 
tion exactly and (ii) scale to large systems needed to treat 
systems near second-order phase transitions. For these 
reasons we use quantum Monte Carlo@(QMC) which in 
principle is exact since there is no "fermion sign prob- 
lem" in this particular ID model (so long as there is 
non-zero overlap between our trial function and the true 
gromid state). To our knowledge this is the first QMC 
study of polarization and localization in any system, and 
the first study of the ionic Hubbard model with systems 
large enough to determine quantitatively the nature of 
the transitions and whether or not there exists the spon- 
taneously bond-ordered phase proposed by Fabrizio, et 
al.H. Furthermore, if there are indeed quantum phase 
transitions in the ionic model - whereas it is known that 
there are none in the usual non-ionic Hubbard model - 
then it follows that one must address the issue: Is a crit- 
ical degree of ionicity required, or is the usual Hubbard 
model unstable to infinitesimal ionic perturbations? It is 
known0~0 that the usual Hubbard model is unstable to 
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dimerization at all U. Thus a second question is: does 
this instability play a fundamental role in stabilizing the 
bond-ordered state? 

The organization of the paper is as follows. In Sec. II, 
we introduce the model studied in this paper. In certain 
cases, depending upon the parameters of the Hamilto- 
nian, this model is exactly soluble. We discuss the rele- 
vance of these solutions to the more general case studied 
in this paper. In Sec. Ill formulas for evaluating the 
electronic polarization and localization are presented. In 
Sec IV, we introduce the Quantum Monte Carlo (QMC) 
methods employed to evaluate expectation values and we 
describe their respective limitations. These are Varia- 
tional and Green's Function Monte Carlo algorithms and 
the "forward walking" method for computing expecta- 
tion values of operators that do not commute with the 
Hamiltonian. Our results are presented in Sec. V and 
comparisons are made with previous studies using ex- 
act diagonalization and Hartree-Fock. In Section VI, we 
discuss the differences between our results and previous 
studies and the consequences of our new findings. 

II. THE MODEL 

The generalized ionic Hubbard Hamiltonianil is de- 
fined by 

H = Hoito, U) + i?Ion(A) + IlDim(x), (1) 

where Ho is the Hamiltonian of the usual Hubbard model 

L 

Z,fT i — 1 

Here c^^ „{ci^a) creates (destroys) an electron of spin a 
on site s while fii^c! = c[ cr'^i,a is the density operator of 
electrons of spin a on site i. This system is an idealized 
model of a chain of atoms that can have at most 2 elec- 
trons of opposite spin per atom. The magnitude of the 
matrix element [to) controls the strength of covalency in 



the centrosymmetric lattice and determines the width of 
the energy band in the non-interacting limit. Interactions 
are included only for electrons that occupy the same site, 
and the strength of electronic correlation is determined 
by the ratio of U /to- 
The ionic term, 

i7ion(A) = A^(~1)X^, (3) 

consists of an on-site energy(± A) that alternates between 
neighboring sites, which is intended to model the electro- 
static potential of cations and anions in an ionic material. 
By adjusting the ratio of <o/A, we can vary the degree 
of covalency and ionicity to levels similar to those of real 
insulating systems. 

Although we will not study dimerization, per se, it is 
crucial to include a dimer term that breaks the inversion 
symmetry and is defined with the Su-Schrieffer-Heeger 
forixii 

HT)im{x)=dB. (4) 

Here S — ax denotes a dimerization term in the Hamilto- 
nian {ti = <o(l + (— 1)'(5)) that incorporates the effect of 
alternately displacing the atoms ztx from their equilib- 
rium positions {R{i)o = ia) and a is the linear electron 
phonon coupling constant. The operator B is the "bond 
order" operator 

i 

= I E(-i)^ [E(4i,.c^,'^ + 4,c,+i..)], (5) 

i a 

which is a staggered hopping operator, the expectation 
value of which is the average difference in kinetic energy 
associated with the two bonds in a unit cell. Here N is 
the number of sites, iV/2, the number of cells, and Bi the 
strength of the z*'* bond. (Fabrizio, et al., refer to this 
as a "dimerization" operator; however, we will use the 
term "bond order" @, since it denotes a property of the 
electronic state and may occur even if the lattice is not 
dimerized.) 
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Exact analytic solutions for Eq. ^ exist in several lim- 
iting cases. In the non-interacting case {U = 0), the 
electrons fill the lowest energy band (E(k)) 



E{k) = ±[A'^ + Atlcos^ik) + A{axy sin^ik)]^/^ 



(6) 



up to the Fermi k- vector ±-^. In the case A = a; = there 
is no gap at the Fermi surface and the system is metallic, 
but for any finite A or x a gap is opened at the Fermi 
surface and the system is a band insulator. If A = 
and we perturb the system by adjusting x ^ the lattice 
is known to suffer the famous Peierls instabilityElS and 
energetically favors dimerization. 

Exact solutions in the presence of correlation (U ^ 0) 
are restricted to cases in which (i) there is no intrasite 
coupling (t = 0); (ii) there is a large displacement such 
that 6 = 1 and the lattice is completely deformed into 
an array of independent dimers; or (iii) the case of the 
usual Hubbard model where there is no ionic potential 
or lattice deformation (A = (5 = 0) for which there are 
exact analytic solutions for all [70. In the last case, 
the exact solution predicts that at half-filling the system 
becomes a Mott insulator for any non-zero U. There is 
no change of symmetry from the case of [/ = (which 
is a metal) and at "very large" U/tg the system reduces 
to the Heisenberg spin model, which also has no long 
range order or spin gap in one dimension. For large U 
the exchange coupling of the mapped spin model is J = 
4i^C//(C/^ -4A2). The MI and BI regimes are commonly 
distinguished from one another in literature on the basis 
of spin-charge separationil. In both cases there is a gap 
to charge excitations but in the MI state the spin gap is 
zero while in the BI state both spin and charge gaps are 
non-zero. 

The limiting cases (i) and (ii) are also instructive for 
our purposes. In the former {to = 0) there is a transition 
at A = J7 from a singlet state with two electrons on 
the site with on site energy —A, which is like a band 
insulator, to a state with one electron per site which has 



a spin on each site and is like a Mott insulator. Thus one 
might expect a transition from the BI state to some other 
phase as U is increased even if <o 7^ 0. The second case 
(ii) with (5 = 1 and to ^ always leads to a singlet ground 
state for the isolated dimers, which relates to the known 
result that one has a singlet state with a gap for both spin 
and charge excitations for any degree of dimerization. 
Thus one can ask; does a transition occur from the BI to 
MI regime as (5 — > for [/ ^ 0? Is there a spontaneousil 
bond-ordered phase? We shall test these ideas with our 
QMC simulations applied to the general case where there 
are no exact analytic solutions. 

III. ELECTRONIC POLARIZATION AND 
LOCALIZATION 

The issues associated with calculating the electric po- 
larization in an extended system have a long, torturous 
historyilH. Only recently have formulas been devised 
that express the polarization and localization of electrons 
directly in terms of the ground state wavefunctiorii'0@. 
One type of formulation measures the change in polar- 
ization as a Berry's phase obtained by integrating over 
twisted boundary conditions and an adiabatic parameter 
that characterizes the evolution of the system as it moves 
from one state to anotheil§'0. This approach has also 
been extended to localization in an independent particle 
formulatioii0 and recently in a many-body formalismEl. 
An alternative approach has been developed by Resta 
and Sorella00 and otherfSQ, who expressed the elec- 
tronic polarization and localization in terms of the ex- 
pectation value of a complex operator 



< Z >=< e 



>=< I I e^^^''^ >, 



(7) 



where the average is taken with respect to a truly corre- 
lated many body wavefunction utilizing periodic bound- 
ary conditions (PBC) sampled using one of the quantum 
Monte Carlo techniques discussed in section IV. In terms 
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of < Z > the polarization of the many body ground state 
can be expressed as 

< APei >= hm ^ Im In <Z>, (8) 

and a measure of the electronic delocalization is given by 

L 



< AX^ >= lim -( — )^ln| <Z > 



(9) 



These expressions are exact only in the limit of an in- 
finitely large system, and in practice one measures each 
for increasingly larger supercells until convergence is met. 
Recently Souza et ali have shown that Eqs || and ^ are in 
fact valid in a correlated many-body system and related 
this formulation to that using twisted boundary condi- 
tions. They also demonstrated that the formulas relate 
directly to measurable fluctuations of the polarization, 
thus validating the two formulas as direct measures of 
electronic polarization and delocalization. 

To our knowledge the present work is the first study 
of polarization and localization on large systems with 
fully-correlated many-body wavefunctions sampled using 
QMC. Previous work has been limited to exact diagonal- 
ization studies on small systems or mean field methods 
such as DFT and HF. For our work we use quantum 
Monte Carlo methods techniques with Eqs || and |9| be- 
cause these are directly in the form of expectation values 
of quantities using wavefunctions that have the usual pe- 
riodic boundary conditions. This is a great advantage in 
QMC since we can use the same methods developed for 
other problemsH The approach using twisted bound- 
ary conditions would require a change in the algorithms, 
in particular the adoption of a "fixed-phase"@jll rather 
than a fixed node method. Such an approach would have 
important advantages, the most significant that it would 
allow calculations of polarization and localization to be 
done on smaller supercellS. There are other reasons 
that we prefer to use the standard boundary conditions: 
we shall see that very large cells are readily handled in 
QMC and furthermore the ability to work with large sys- 
tems is very important in conclusions on the nature of the 



phase transitions in this study. 



IV. QUANTUM MONTE CARLO 



Quantum Monte Carlo (QMC) methoddfJ'B make it 
possible to evaluate expectation values of operators in 
many-body systems by stochastically sampling a proba- 
bility distribution. In this paper we focus on two meth- 
ods. Variational Monte Carlo (VMC) and Greens Func- 
tion Monte Carlo (GFMC), that can be used to deter- 
mine properties at temperature equal zero. The space 
of integration (R) is the set of all the electronic coordi- 
nates {r*!, ... , rw}, which is sampled by "walkers" which 
denote a set of configurations R. A random walk is gen- 
erated by starting from an initial configuration Rq, from 
which new configurations are generated by successively 
stepping to new random configurations, e.g., using a gen- 
eralized Metropolis methodEl. This is done by accepting 
or rejecting new configurations at each step based upon 
a chosen acceptance function (P(R)). After a period 
of time the walk will stabilize such that the set of con- 
figurations visited {R} will be distributed according to 
P(R). VMC measures expectation values by uniformly 
averaging over the configurations visited by the Metropo- 
lis algorithm where as in GFMC the average is weighted 
according to R. 

A. Variational Monte Carlo 

VMC measures expectation values of a variational trial 
wavefunction ( 5't({q!},R) ), where {a} denotes a set 
of parameters that can be optimized. Averages for an 
arbitrary operator O are obtained by sampling 

/ *T({a},R)6*T({a},R)dR 



VMC 



/|*T(W,R)|2ciR 
/|vI/T({a},R)p6L(R)dR 



(10) 



/|M'T(W,R)PdR 

the local form of Ol, defined as 0^t(R)/^'t(R.), over a 
set of points ({R}) distributed according to the modulus 
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of the wavefunction. The {R} are obtained by choos- 
ing I^'tP as the acceptance function in a gencraUzed 
MetropoUs algorithm. VMC is easy to implement but 
is limited in accuracy by the form of the adopted wave- 
function. In our work has the Gutzwiller formS 



Trial 



Jastrow Term (ff) 



which is a product of Slater determinants for each 
spin (thus guaranteeing that the wavefunction is anti- 
symmetric) and a two body Jastrow correlation function 
that reduces the amplitude of configurations with doubly 
occupied sites for < g < 1, thus lowering the interaction 
energy. The single body portion of Eq. |l| is parameterized 
by A' and S' , which means the orbitals used to construct 
the Slater determinants are obtained by diagonalizing 
the non-interacting {U = 0) portion of the Hamiltonian 
{H{A' ,d')) and adjusting {A',S') to optimal values that 
minimize the energy in Eq. wrt. 5't(5, A' , S'). 

B. Green's Function Monte Carlo (GFMC) for 
Discrete Systems 

GFMC starts with the optimized VMC wavefunction 
^'t(5, A', 6') upon which a projection is applied to obtain 
an improved ground state. To illustrate the principles 
upon which this method depends, one can expand '^t in 
terms of the eigenstates of H. Then the imaginary time 
propagator acting upon 'i'x has the form 

^ — ^ r — >co 

n 

Note that the exact ground state, ^'o, can only be ob- 
tained so long as it has non-zero overlap with '^t- 

The following is a summary of the method developed 
by Haaf et alEll some of which is used in the next section. 
For lattices this projection scheme takes advantage of 
the fact the spectrum of H is bound such that one can 
use a Green's function projection with no finite-time-step 



erro 



^-t{H-Eo) 



[1 - At{H - Eo)] 



NAt 



(12) 



The propagator acting upon the trial wavefunction now 
becomes 



I* 



N\ 



= [l-AriH-EoT I^t) 



By inserting the identity operator in the real space con- 
figuration basis (R = {ri.t, • • • , ^n,i}) 

R 

between successive applications of the projection opera- 
tor and multiplying both sides by (Rn| 

*^(Rn) = 

E (Rn|[1- Ar(H-Eo )]|Rn-i) 

R]sj_i,... ,Ro 

(Rn-i|...|Ro) (RoI^t) 



(13) 



we obtain an expression for the wavefunction after N 
steps in imaginary time. If the time step At is sufficiently 
small (R| [1- At(H-Eo)] |R') > and can be interpreted 
as a probability. Using this probabilistic interpretation, 



the sum in Eq. |13| above is evaluated using Metropolis. 
Multiplying and dividing by (Rj^* Trial), Eq. ^ above can 
be importance samplec0 as 



*~(R 



N 



N 



E *?'(Rn) nG(R.i,Ri-i)*T(Ro), (14) 

where 

G(Ri, Ri_i) = (Ri|[l - Ar(H - Eo)]|Ri-i) 



*T(Ri-l) 



(15) 



Since the G(Ri,Ri_i) are not normalized to one, they 
can not be interpreted directly as a probability. This is 
remedied by expressing G'(Ri,Ri_i) as 



G(Ri,Ri_i) = m(Ri,Ri_i) p(Ri,Ri_i), 



(16) 
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where p(Ri, Ri-i) is identified as the probabihty of mov- 
ing from Rj to Ri_i and given by 



p(Ri,Ri_i) = |G(Ri,Ri_i)l/|m(Ri,Ri_i)|, 



(17) 



while the weight TO(Ri,Ri_i) normahzes p such that 
X;i{,p(Ri,Ri-i) = 1 and is 

m(Ri,Ri_i) = sign(G(Ri,Ri_i) )^ |G(R',Ri_i)|. 

R' 

The nodal structure of the ground state divides config- 
uration space into regions in which 5*71 (R) is positive or 
negative, so that G(Ri,Ri_i) changes sign upon cross- 
ing the nodal surface in configuration space. Crossing 
the nodal surface by a walker causes difficulties in Monte 
Carlo sampling since the weight of a walker must be pos- 
itive definite if it is to be interpreted in a probabilistic 
manner. In general, one must make some approximation 
to remedy this problem, by fixing the sign of G in the 
Monte Carlo sampling; this is referred to as the "fixed 
node approximation" , which has been described for lat- 
tice problems by ten Haff, et akEl 

In the generalized Hubbard model considered here, 
QMC is exact because: (i) the only nodes of the ground 
state wavefunction are the points where two electrons of 
the same spin cross, (ii) the nodes are the same as in 
the trial function which automatically obeys this condi- 
tion, and (iii) the Monte Carlo sampling is restricted to 
a nodal region in which the sign of G(Ri, Ri-i) is fixed. 
The last condition is realized in the present work because 
each move involves only one electron moving one site at a 
time; we never reach the nodal surface since neighboring 
R in which a site is doubly occupied by two electrons of 
the same spin are not allowed. Thus our algorithm sam- 
ples one nodal region (either positive or negative) which 
is sufficient, since they are identical due to the antisym- 
metry of the wavefunction. 

Implementation of the above method is as follows. A 
VMC calculation is performed which supplies a number 
of walkers {R} initially distributed according to |^tP- 



Each of these are then randomly walked along a path in 
configuration space using p(R, R') as the Metropolis ac- 
ceptance function of moving from R to R'. Each step is 
weighted by m(R, R') such that the i*'* walker's accumu- 
lated weight is 

TV 



..N 



Expectation values for an arbitrary operator O after N 
projections of the green's function are measured by aver- 
aging the weighted local form of O of each walker 

(^-tIOI*^) _ E,OL(RN)wf' 



GFMC 



TV\ 



(18) 



Averages in GFMC equal the ground state expecta- 
tion value only for those operators which commute with 



H because the inner product Eq. 18 is a "Mixed Estima- 



tor" between (4't| and |^o}- Operators that commute 
with H share the same eigenstates and the operator in 
Eq. |l8| can be considered to act to right on ^'o, thus 
returning the ground state and cancelling the normaliza- 
tion of the denominator. Conversely operators that do 
not commute with H have different eigenstates and thus 
do not cancel the normalization of the denominator in 
Eq. |l^. Consequently GFMC does not produce exact re- 
sults for these operators; such expectation values will be 
addressed later. 

C. Test of GFMC on Ordinary Hubbard Model 

The accuracy with which the energy can be measured 
in GFMC and the magnitude of finite size effects can be 
addressed by comparing with exact results for the usual 
Hubbard Model at 1/2 filling, which have been evaluated 
by Hashimoto0 for finite systems of AN + 2 sites using 
periodic boundary conditions. In Fig |^ the differences 
in energy between lattices of size L and the thermody- 
namic limit is plotted for two cases U/t — 1.25, 5.0. The 
finite size effects at typical U « 2.4 are of order 0.0001 1 
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for a supercell of 82 sites. Thus we do not anticipate 
any difficulty in calculating the energy except in cases 
where there is a much longer correlation length than in 
the usual Hubbard Model, e.g., near a phase transition 
where correlation lengths diverge. 
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FIG. 1. In the lower figure E{AN+2)-E{oo) is plotted where 

A'' = 2,3,... ,11 and infinite system estimates are those of Lieb 
and Wu. The lines are exact results from HashimotoS and the 
symbols are the QMC estimates for i) U — 5 (diamonds) and 
ii) U = 1.25 (squares). In the top figure the energy difference 
between the QMC and exact results is plotted vs L. 



D. Stability of the Ordinary Hubbard Model (A = 0) 

For comparison to our work later it is useful also to 
study the dimerizcd Hubbard model with S ^ 0. Work 
on related issues in the past two decades has verified 
early theoretical predictions0 that electron correlation 
enhances the Peierls instability of the non-interacting 
Hubbard model as 5 ^ 0. Using the Hellman-Feynman 
theorem the bond order {B) can be identified as the first 
derivative of the energy wrt the lattice distortion ax or 
S. The bond order susceptibility or the second deriva- 



tive of the energy wrt. 5 has a logarithmic divergence as 
which is referred to as the Peierls instability. 
The energy near 6 = varies asii 



E{6 = 0) + AS''/\n{6), 



(19) 



where the amplitude A and 7 are dependent upon the 
strength of electron correlation. For U = A is pro- 
portional to to and 7 = 2, and for U/to << 1 varia- 
tional methods suggest the same results. In the strongly 
correlated regime the lattice can be mapped onto a ID 
Heisenberg lattice where A is proportional to 4t^/t/ and 
7 = 4/3. Although the instability is enhanced at large U, 
the effect is more difficult to observe since the electronic 
energy is much smaller. 

In our studies we consider small ionic deviations {S ^ 
0) from the usual Hubbard model for U = 2.4. The QMC 
energy and bond order are plotted vs 6 in Fig ^ for an 82 
site lattice. The GFMC energy was fit to Eq 19 using a 
non-linear least squares routine. The parameters of the 
fit are E{S = 0) = -0.777589(24), A = 1.48(17), and 
7 = 1.29(3) and give a reduced chi square of 1.58. This 
data agrees quite will with that of Black and Emerj^H 
who observed 7 = 4/3 in the ID Heisenberg model. The 
energy of the symmetric lattice is within error bars of 
the exact thermodynamic limit of —0.77762. The diver- 
gence of the lattice's susceptibility of the lattice to bond 
ordering can be observed in Fig |2[ as the level of dis- 
tortion approaches zero the bond order approaches the 
origin with infinite slope. 
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FIG. 2. Ground state energy and bond order vs lattice dis- 
tortion 5 for U = 2.4 in the usual Hubbard model. The energy 
was fit to the function E{S = 0) + AS'' / ln{5) using a non-linear 
least squares routine. 



E. Expectation Values and Forward Walking 

As noted before, GFMC does not produce exact expec- 
tation values for operators that do not commute with H. 
There are several ways to improve upon the GFMC mixed 
estimator for such expectation values. One is an approx- 
imation that is valid so long as the VMC and GFMC 
averages are close to one another. Expressing I^I^o} as 
+ 1'^^'} and taking the inner product, the ground 
state expectation value can be expressed afl 

(^'0|O|*0) ~ 2 {d)GFMC - {d)vMC + 0{d^^). 

(20) 

However, this approximation breaks down whenever the 
VMC trial wavefunction is not a good approximation to 

The exact ground state expectation value of any oper- 
ator (O) can be found if the mixed expression Eq. ^ is 
replaced by one involving the exact wavefunction in both 



the bra and ket 

(*t|[1 - AriH - Eo)]''6[l - Ar(i? - Eo)]^\^t) 



{^t\[1 - At{H - Eo)]M[l - At{H - Eo)]N\^t) 



(21) 



This can be accomplished by "forward walking"ii, which 
can be simply expressed in terms of the GFMC method 
previously discussed. The same methods and terminol- 
ogy used in GFMC are also applicable here. Inserting 
the identity operator between each projection and using 
importance sampling Eq. ^ can be rewritten as 

N+M-l 

^ [ n G(Ri+i,Ri)] Ol 

R-N + ivi,--- i—N 

N-1 



[[] G(Ri+i,Ri)]^'^(Ro). (22) 



1=1 

The G(R, R') are sampled as before in terms of a prob- 
ability function (P(R,R')) and weight (A/(R,R')). A 
series of i walkers, initially distributed according to the 
VMC trial function, are stepped along paths ({Ri}) 
in configuration space by Metropolis sampling. After 
N projections the accumulated weight of each {Ri} is 
the product of all steps weights, as defined in Eq. 
The walkers weights are distributed according to the 
mixed probability distribution \E't(Rn)^o(Rn)- The lo- 
cal form of O (Oi(RN) is measured for each walker but 
not averaged as it is in GFMC. The walkers are moved an 
additional M steps in imaginary time over which they ac- 
cumulate post measurement weights {wf^). Averages are 
computed using each walkers accumulated weight before 
and after measuring Oi(RN) 

E.^f[<o.(RN)] 

Although this method is in principle exact, assuming the 
nodal structure of is known, it also has its disadvan- 
tages. In particular, the width of the post-measurement 
weight distribution grows with the number of steps M, 
thereby increasing the fluctuation of the forward walking 
estimates. To achieve a desired level of accuracy addi- 
tional measurements are needed but the error in QMC 
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decreases inversely with the square root of their number. 
Consequently, to obtain the same error as that in GFMC, 
forward walking may require many times more estimates 
in Eq. |2^. This is the limiting factor in applying forward 
walking. 

To illustrate the practicality and usefulness of this 
method, we show in Fig ^ the bond order < S > as 
a function of forward walking for the ccntrosymmetric 
lattice at C/ = 1.8, A = 4/7 and L = 62 sites. We have 
chosen a poor trial wavefunction biased towards bond 
ordering by defining the determinant part of the wave- 
function Eq. |l^ using a Hamiltonian with 5' = 2/35. 
At this particular U the system is a band insulator (as 
shown below), consequently the bond order must be zero; 
however, the average bond order in VMC and GFMC is 
non-zero as a result of using this trial wavefunction. The 
VMC expectation value of the bond order is substantially 
different from and reflects the poor quality of the trail 
state; whereas the GFMC average is closer to the exact 
result but remains far from satisfactory. The results in 
Fig H illustrate the strengths and weaknesses of forward 
walking. The key point is that there is a competition be- 
tween the improvement of the estimate and the growth of 
the statistical errors with projection time. As shown in 
the figure, the method vastly improves the results even 
for very poor wavefunctions without the proper symme- 
try. In general, we use much better trial wavefunctions, 
and so the convergence to the exact result in our work 
below is more rapid than that depicted in Fig ^. 
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FIG. 3. Illustration of forward walking for U — 1.8 and 
{A, S} = {4/7, 0}. The results are those obtained using a wave- 
function of the form in Eq ^ with variational parameters {A, 5} 
— {0.38, 2/35}, which is much worse than a typical optimized 
trial wavefunction used in the present work. 



V. RESULTS FOR IONIC HUBBARD MODEL 

The unit cell for the ionic Hubbard model is composed 
of 2 sites and the Hamiltonian is given in Eq. |l|. In order 
to understand the meaning of the polarization and bond 
order in this system, it is helpful to consider first the 
non-interacting case with U = 0, where one can visualize 
the electronic properties in terms of Wannier functions. 
At zero dimerization (S ~ 0) the Wannier functions are 
centered on the sites whose onsite energy is shifted by -l-A 
and — A. The 2 electrons of opposite spin in each unit cell 
both occupy the lowest energy Wannier function centered 
on the lower energy site. In the dimerized lattice {8 7^ 0), 
the centers are displaced from the sites. The magnitude 
of (5 dictates the amount by which they are off center from 
the lowest energy site while sign of 5 determines whether 
the center is to the left or right of this site. Strictly 
speaking the existence of Wannier functions in correlated 
systems ([/ ^ 0) has yet to be proven but we will continue 
to use the concept of the center of the localized states for 
illustrative purposes. As U increases the electrons find it 
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energetically undesirable to occupy the same site, and in 
the dimerized state the center of the distribution shifts 
further away from the low energy site. The limit of this 
displacement (i.e., polarization) is P = 1/2 since one 
would never favor having more than one electron on the 
higher energy site. Similarly, the limit of the bond order 
is {13) — 2 corresponding to isolated dimers. 

As (5 — ^ there arc only three possibilities. If there is a 
spontaneous breaking of the inversion symmetry the po- 
larization can assume any fractional value between and 
1/2. If there is no breaking of symmetry, there are still 
two possibilities since there are two centers of symmetry: 
the polarization can be or 1/2. If the reference point 
defined to be zero is the usual band insulator where both 
electrons occupy the Wannier function centered on the 
lower energy site, it has been proposedElil'0 that P = 
1/2 corresponds to a Mott insulator with no long range 
order. 

We first report results of our study of the ionic Hub- 
bard model with parameters fixed at the values used 
in previous work0S0, so that direct comparisons can 
be made. The energy scale is set by defining to — 1, 
A/to = 0.5714 and aa/to = 40/7. The previous con- 
clusions with which we will compare are based upon 
exact diagonalization of the many-body Hamiltonian 
in small supercellS'll and Hartree-Fock calculations!!. 
The studyEl using exact diagonalization of 8 site lattices 
with twisted boundary conditions found a jump of 1/2 
in the electronic polarization for 6 — 0, i.e. an electron 
in each unit cell being transported 1/2 lattice constant, 
at a critical value of U (Uc — 2.26). This was inter- 
preted as a transition between BI and MI phases, which 
was supported by Hartree Fock (HF) calculations that 
showed similar behavior at Uc — 2.46. Extrapolations 
using larger cells of 12 siteJ3 find Uc = 2.86, presum- 
ably a more converged value. The key points are: (i) the 
transition point Uc is found to be a metallic point with 
divergent dclocalization; (ii) effective charges diverge and 
change sign at the transition; and (iii) there is no sign 



of the bond-ordered state predicted by Fabrizio, et al. 
This new state would have long range order and break 
the inversion symmetry of the lattice, thus allowing the 
polarization to take any fractional value. 

The present work is based upon the QMC algorithms 
described earlier and the formulas for polarization and 
localization in section III. The first step in applying the 
QMC methods is to find a trial wavefunction that has 
as much overlap with the true ground state as possible. 
This is achieved by optimizing the parameters {g, A' ,S'} 
to minimize the energy. To determine the optimal value 
of g we have used a newly devised technique that sig- 
nificantly reduces the amount of computational effort 
requiredii. Using the optimal Gutzwiller parameter the 
energy of '^rig, A' ,6') for different A' and S' is sampled 
using VMC. We adjust A' and 5' to lower the VMC en- 
ergy and measure it at several points in the neighborhood 
of its minimum. A curve fit is then performed using these 
points to determine the optimal A' and 6' . 

A. Comparison with Exact Diagonalization and 
Hartree-Fock 

Previous studies distorted the lattice by varying de- 
grees and these results provide a basis of comparison 
with QMC. We have measured the polarization of the 
ionic lattice for large {S — 0.08) and small {S — 0.02) lat- 
tice distortions and plotted these with the corresponding 
results of previous studies in Fig ^. Size effects are ac- 
counted for by extrapolating to the thermodynamic limit 
in 1/L; this will be outlined more clearly in the follow- 
ing section. The Lanczos results agree well with those 
of QMC for S — 0.08 considering the fact they were ob- 
tained using 8 site supercells with twisted boundary con- 
ditions. This is in agreement with previous studies using 
exact diagonalizationilS on the usual Hubbard model 
which found that small cells of this size were sufficient to 
reach thermodynamic convergence in the 0.05 < S < 0.1 
regime, whereas convergence with cell size is worse for 



11 



smaller S. (The jump in the polarization found i 
calculations for non-zero 6 is unphysical and arisi 
cause the mean field approximation leads to an 
ferromagnetic ground state. In ID this is strictb 
hibited because quantum fluctuations are strong e: 
to destroy long range order in any continuous qua: 
As the magnitude of the distortion (6) approaches 
difference between QMC and Lanczos becomes gr 
Studies using exact diagonalization and HF observe^ 
the polarization as a function of S tended to be 
critical U and 1/2 above it; consequently the dy: 
charge was observed to change sign upon crossin 
critical point. In the lower plot of Fig ||this is exh 
as a crossing of the curves for P{6). On the couuiaiy 
we observe that the dynamic charge remains the same 
sign for the entire range of U studied (except that the 
sign of dP/ d6 is difficult to establish for large U where 
it is near zero). This difference is attributed to the fact 
that for small 6 = and for U near Uc the electrons 
are very delocalized and correlation lengths exceed the 
cell siz E. Integrating over twisted boundary conditions 
provides thermodynamically quenched expectation val- 
ues so long as the Wannier functions of the eigenstates 
^'fe have vanishing overlap. Near the critical point the 
^'fe, obtained by exact diagonalization of 8 and 12 site 
rings, do overlap significantly; consequently, regardless 
of the number of k-space points averaged over by Resta 
and Sorella, the polarization will not converge to that of 
QMC. 




U/to 

FIG. 4. Exact QMC measurements of polarization ( upper 
figure ) in comparison with previous Lanczos and HF results 
( lower figure ). Results are illustrated for staggered transfer 
integrals of 1± 0.02 ( squares ) and 1±0.08 ( circles ). The HF 
results are depicted by the dashed line in the lower figure. 

The average polarization in Fig ^ was measured using 
forward walking. For the polarization at these values of 
6, forward walking is not essential and the same results 
can be obtained using Eq ^ However, localization is 
more sensitive to inaccuracies in the wavefunction and 
accurate expectation values are provided only by using 
forward projection even at these relatively large values 
of 5. As the level of dimerization is reduced such that 
6—^0 the extrapolation technique breaks down and only 
forward projection can provide accurate estimates for po- 
larization and localization. Consequently we only report 
in this paper those QMC results obtained using forward 
projection. 

The observation of a topological transition in the work 
of Resta and Sorella, and Guidopoulos, et al., is based 
upon their finding that as (5 — > the polarization jumps 
discontinuously from to 1/2 at a critical Uc = 2.26 (or 
Uc = 2.86). This means that the dynamic charge (Z), 
defined as dP/dx\s=o, diverges and changes sign at Uc as 
6^0. In latter wor ki Resta and Sorella showed that 
Uc is a metallic point where the electronic localization 
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length (< A^X >) diverges, and Guidopoulos, et al., 
found energy gaps that extrapolated to zero. We will 
compare these results with our work below. 

B. Phase transition to Bond-ordered State 

We have measured the forward walking estimators for 
P, < B > and < A'^X > and taken the limit of (5 ^ 
to study the nature of the quantum phase transition. 
The formulas used to obtain expectation values for po- 
larization and localization are only accurate in the limit 
L oo. This limit is taken by fitting measurements at 
finite L to a linear least squares fit in (l/L)'^ and ex- 
trapolating to 0. We have found 7 = 1 to accurately 
account for the finite size effects of P and 7 = 2 for 
< A'^X >. This scaling has only been found appropriate 
upon increasing the supercell size above a critical thresh- 
old which depends on the proximity of the metallic state. 
The accuracy of the finite size corrections to P are il- 
lustrated in Fig H at [/ = 2.7 for different magnitudes 
of S. The data in Fig | was collected near the critical 
point of the phase transition, where size effects are large 
and must be treated accurately. If the system is suffi- 
ciently far from such a critical point, size effects are k 
pronounced and there is a more rapid convergence to t 
thermodynamic limit. 




FIG. 5. The points represent the QMC P ob- 
tained using forward walking for system sizes of 
80, 100, 140 and 200 sites for different levefs of dimer- 
ization S = {0.0028,0.0056,0.0085,0.0114} in ascend- 
ing order. 

Using the infinite L estimates for the polariza- 
tion and localization on lattices dimerized hy 6 = 
{0.0028,0.0056,0.0085,0.0114} we have performed a lin- 
ear least squares fit and extrapolated to the centrosym- 
metric limit {6 — 0). This method makes the assump- 
tion that the response of the lattice to dimerization is 
linear. However, near the phase transition non-linearity 
will cause this to break down. In Fig |^ we have plotted 
the polarization and localization of the ionic model for 
different magnitudes of dimerization. 




FIG. 6. P{L = 00) and < A^X > for various 
levels of dimerization (5). The extrapolated cen- 
tro-symmetric polarization and localization is repre- 
sented by the points with error bars. 
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The phase transition we observe using QMC differs 
from the topological transition found using exact diago- 
nalization and HF. As previously mentioned, in the BI 
and MI phases the polarization is restricted to or 1/2. 
Resta and Sorella identified the shift from to 1/2 as 
the signature of a BI ^ MI transition. However, we 
observe that P takes a continuous range of values in the 
centrosymmetric limit which cannot occur in either the 
BI or the MI phase. This can only occur if the global in- 
version symmetry of the lattice is broken by a long range 
bond ordered state, predicted by Fabrizio et aS on the 
basis of field theory arguments in which he mapped the 
Hamiltonian onto two Ising spin models. The order pa- 
rameter of this phase transition is the average bond or- 
der function < _B >, where B is given in Eq|. The 
spontaneous bond order of the centrosymmetric lattice is 
obtained by extrapolating to (5 = the {B) of the same 
distorted lattices as before. We fixed the supercell size to 
142 sites and found the consequent size effects are within 
order of the error with which we can measure the bond 
order. These results are plotted in Fig 0. 



fflffl 



£8 



m 
m 



eeeee ' <° 



U/to 

FIG. 7. QMC Bond Order of centro-symmetric lat- 
tice for A/to = 0.5714 and L = 142. Results obtained 
by extrapolating bond order on distorted lattices of 
5 = {0.0028, 0.0056, 0.0085, 0.0114} to 5 = 0. 

We have attempted to classify the quantum phase tran- 
sition by fitting the polarization and bond order of the 



centrosymmetric lattice to a function of the form 

A[U~Uc]^, (24) 

where ^ is the critical exponent and determines the 
universality class of the transition. A non-linear least 
squares routine was used to fit the data, with fitted pa- 
rameters Uc, A, and ^ listed in Table I. In Fig ^ the data 
for P and < B > and the corresponding fits are plotted. 
Both quantities behave similarly near the critical point 
and the Uc of each is nearly identical. We find for P 
and < B > are near 1/2, the expected mean field ex- 
ponent. On the other hand, Fabrizio, et al.E!, predicted 
that the transition is of the Ising universality class and 
thus ^ should be 1/8. We do not know whether the dif- 
ference is real or it is simply due to the possibility that 
the range oi U — Uc over which the scaling belongs to 
the universality class is too small for us to observe in the 
present work. 
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FIG. 8. QMC polarization and bond order near the 
critical point and their relative fits to Eq 

Alternatively, the existence of the bond ordered state 
can be observed by directly studying the symmetric lat- 
tice without any lattice distortion. Quantum Phase tran- 
sitions (QPT) are characterized by a symmetry breaking 
that occurs in the thermodynamic limit. Below Uc the 
lattice is a band insulator with no bond order but above 
Uc the electrons will spontaneously choose to bond order 
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with ± I (B) I . There are two such states characterized by 
the same magnitude but opposite sign of the bond order. 
For any finite system the ground state remains a hnear 
combination of both. However in the hmit L —>■ oo one 
of these is arbitrarily chosen as the ground state. Even 
though the QMC simulations of the symmetric lattice 
for U > Uc measure zero bond order and polarization for 
long simulations, the imaginary time evolution of the sim- 
ulations clearly depict the projected ground state mov- 
ing from one of these bond ordered states to the other. 
This phase separation gives rise to large auto correlation 
times. The evolution of the bond order and polarization 
in imaginary time are illustrated in Fig ^ for U — 3.45 
and L — 60 sites. As the ground state moves between 
BO states of opposite symmetry both the polarization 
and dimerization are observed to change sign. This pro- 
vides an alternative method of detecting the existence of 
the BO state. 
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FIG. 9. Illustration of the phase separation in the 

symmetric case {S = 0,A = 0.5714,[/ = 3.45) of the 
polarization and bond order. The lines depict the mea- 
surements over which averages are obtained in QMC. 



At large U one might expect that the ionic ID Hub- 
bard model is a Mott Insulator. Fabrizio, et al.El, pre- 
dict the existence of a Kosterlitz-Thouless transition for 
large U/ A at which the lattice becomes a Mott Insulator. 



Such a transition would be characterized by polarization 
of 1/2 and no bond order. On the contrary we do not ob- 
serve either for any U considered which included values 
up to C/ = 10. The bond order does diminish but ap- 
pears to asymptotically approach and the polarization 
appears to converge to 1/2 only in the limit of [/ — > oo. 
Yet working with such strongly correlated systems has 
the disadvantages that (i) fluctuations of the local esti- 
mators increase due to greater inaccuracies in the trial 
wavefunction and (ii) forward walking works so long as 
the trial wavefunction has some overlap with the exact 
ground state and as U increases overlap with the exact 
ground state diminishes. 

C. Phase transition as a function of ionicity A 

An alternative approach to study the phase transi- 
tion(s) is to diminish the ionic potential A while keeping 
U fixed, so that the ratio U /A increases. We have fixed 
the strength of electron correlation to U = 2.4 and stud- 
ied the bond order and polarization for < A < 0.5714. 
The behavior of the centrosymmetric lattice is inferred 
using two approaches (i) extrapolating results obtained 
on lattices with 5 Q and (ii) looking for evidence of 
phase separation in the symmetric case. Fig |l^ shows the 
results of the first approach for a fixed supercell length 
of 142 sites. In the first we've neglected size effects and 
fixed the super cell length to 142 sites. At large A the 
single body contribution to the Hamiltonian is the dom- 
inant term and the lattice is a band insulator. Conse- 
quently the bond order and polarization are 0. How- 
ever, as A a transition occurs to a BO state where 
the bond order is non-zero and the polarization assumes 
values between and 1/2 as before. (The transition is 
rounded at this fixed cell length.) The bond order in the 
(5—^0 limit is shown by the dotted line in Fig 10. These 



results were obtained by linearly extrapolating the bond 
order at finite S. (No extrapolation was performed for 
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the polarization since it is sensitive to size eilects tha 
were addressed in the previous section.) 

Our results indicate the bond-order state exists at al 
values of A 7^ studied. The finite value of the bon( 



order for (5 — > shown in Fig 10 contrasts sharply wit! 
the vanishing of the bond order as 5 for the non 
ionic Hubbard model (A = 0) as shown in Fig|^. At A = 
6 = 0, we always find < i? >= and polarization equa 
1/2 as they must be for a MI state with no long rang 
order. However, our QMC simulations of the symmetri 
case ((5 = and U = 2.4) at the smallest value of the ioni 
potential studied A = 0.0716 reveal two BO states phaS' 
separating in imaginary time qualitatively the same a 
shown in Fig|9| Thus from our studies, there is no sign o 
a second transition to a MI state as proposed by Fabrizio, 
et al.il, and the long range bond ordered state in Fig |l^ 
appears to exist for any finite A ^ 0. 

This implies that the MI state in ID exists only within 
the usual Hubbard model and in ionic Hubbard lattices 
only in the limit U = oo. At large U the ionic Hubbard 
model has been mapped onto the Heisenberg spin model. 
The present finding suggests that such a mapping may 
be insufficient for Hubbard models with ionic potentials 
and that terms ignored or considered small possibly play 
a fundamental role. 
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FIG. 10. Polarization and Bond Order vs Ionic potential A 
for S = {0.0028,0.0057,0.0114,0.0171} and U = 2.4. The ex- 
trapolated bond order of the centro-symmetric lattice is de- 
noted by the dotted line. 



VI. LONG RANGE ORDER AS INFERRED BY 
BOND ORDER CORRELATION FUNCTION 

Existence of a long range bond ordered state can be 
inferred by measuring the bond order correlation function 
(5s (»'))• We define .gs(r) as 



9B{r) = l{Y. ^»^'+' ) 



(25) 



where Bi is defined in Eq. || and is the strength of the 
i*^ bond of the lattice. If the BO state exists then 
this correlation function would be staggered as a con- 
sequence of the periodic arrangement of dominant and 
weak bonds. In Fig. ^^^(r) is plotted for 4 separate 
cases: (i) A = 0.5714, U = 1.2, L = 60 (u) A = 0.1432, 
U = 2.4, L = 122 (in) A = 0, U = 2.5, L = 122 and 
(iv) A = 0.5714, U = 3.45, L = 60. The first case cor- 
responds to the band insulating regime in which gsi^) 
exponentially approaches a constant, confirming the lack 
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of any long range ordered phase. Conversely in the last 
case, which corresponds to the system in Fig |9| that ex- 
hibited phase separation, its clearly visible that gBif) is 
staggered, signifying the presence of a long range bond 
ordered state. Finite size effects in each of these cases 
were determined to be miniscule and small systems were 
deemed sufficient to measm'e gsif)- The second case 
corresponds to diminishing A so as to move the system 
towards the established Mott State of the usual Hubbard 
model. At this point in the phase diagram the wells of 
the bimodal distribution are weakly defined; thus making 
it extremely difficult to observe the phase separation of 
the bond order parameter directly. In contrast, the bond 
order correlation function is clearly staggered, though to 
a lesser degree than that of the later case. Case iii) is 
the Mott state of the usual Hubbard model. The stag- 
gered behavior of g sir) does not approach a finite limit 
at large r, but rather tends to in a fashion that ap- 
pears to be a power law ; contrary to the exponential 
convergence observed in the BI regime. Comparison of 
cases ii) and iii) shows that in each case the staggered 
behavior of gB (r) is longer ranged than in the band insu- 
lating and strongly bond ordered cases; this exemplifies 
the difficulty of measuring the bond order or the phase 
separation of this generalized model as the ionic potential 
tends to zero. 
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FIG. 11. Bond order correlation functions for i) band insulat- 
ing ii) weakly bond ordered iii) mott insulating and iv) strongly 
bond ordered regimes. See text. 

The staggered nature of gsir) can be used to deter- 
mine the mean bond order of the lattice. Defining two 
new variables 

A5B(r) = [<?s(r)-5B(r + l)].(-ir (26) 



and 



g^{r) = 9B{r)+gB{r + l) ^^^^ 



and substituting Eq ^ in place of gn^f) we can relate 
these two quantities in terms of measurable quantities. 
At large r the i*'* and j*-^ bonds are uncorrelated and 
Agsir) is the RMS bond order of the lattice whereas 
gsir) is the square of the average bond strength. 

In the limit of large r the average bond order (B) can 
be expressed in terms of Ags^r) as: 



V2 Agsir) = {B){r) 



This estimate is exact when r ^ 



{B). 



where 



(28) 



IS 



the correlation length. Fig shows {B){r) plotted vs r 
for the same cases as in Fig |ll|. The strongly bond or- 
dered system converges to an estimate of the bond order 
that is remarkably close to that obtained by extrapo- 
lating from distorted lattices (0.45). The weakly bond 
ordered case appears to converge to a value near 0.18 
which is in reasonable agreement with the extrapolated 
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value of 0.2263(46). The band insulating estimate rapidly 
approaches as a function of r, whereas the Mott insulat- 
ing system appears to approach in a power law fashion. 
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FIG. 12. {B){r) vs r for i) band insulating ii) weakly bond 
ordered iii) mott insulating and iv) strongly bond ordered 
regimes. 



VII. DISCUSSION 

One of the primary results of the present work is the 
quantitative demonstration of the stability of the bond- 
ordered phase for interaction U above a critical value 
Uc{A.) for any non-zero A. Most of our work was carried 
out through dimerizing the lattice by 5 and examining 
the S ^ limit. There were two reasons for this: (i) this 
is an aid in the actual calculations which are stabilized 
by the applied bias, and (ii) the variation with dimeriza- 
tion 5 is important in and of itself. Regarding the second 
point, it is well known that the ordinary non-ionic Hub- 
bard model is unstable to dimerization, with a logarith- 
mic Peierls instability at [/ = that becomes a stronger 
fractional power law instability at large Our 
work shows that as a function of U/A the ionic Hubbard 
model undergoes a phase transition from a stable non- 
dimerized BI phase to a correlated phase in which the 



instability is more severe than in the non-ionic Hubbard 
model. This is evident in comparison of Fig. 2 with Figs. 
5 and 7. In the former for the non-ionic case, the de- 
crease of the bond order with S is clearly observed and 
is consistent with previous theoretical predictions of the 
power law form. However, in all the calculations for the 
ionic model for U /A above the critical value, the average 
bond order < i? > is found to extrapolate to a non-zero 
value. This is observed even for S much smaller than 
previous studies. From this evidence alone there are two 
possibilities: 1) the BO phase with broken symmetry is 
stable at zero dimerization, or 2) there is non-analytic 
behavior as (5 — > which is even stronger than that for 
the non-ionic Hubbard model. 

This result is sufficient to draw conclusions about real 
ID systems in which the sites are allowed to dimerize if 
this leads to lower energy. In cither scenario above dimer- 
ization would always occur (except in the BI phase). In 
the former case the BO phase would occur spontaneously 
and by symmetry there would always be an accompany- 
ing lattice distortion. In the second scenario, dimeriza- 
tion would occur and lead to bond order. The symmetric 
Mott insulator would never occur and the only transition 
would be from a BI phase to a dimerized BO phase. 

At this point we can compare with experiment on ID 
materials. Experimental works by Torrance, et al. ,0 ob- 
served a second order transitions between neutral (BI) 
and ionic (BO) states in organic charge transfer solids. 
The transition occurs upon applying pressure over a wide 
range of temperatures and was attributed to the rise in 
Madelung energy of the crystal. No state synonymous to 
the Mott state was observed. 

In addition, however, calculations on the centrosym- 
metric lattice show directly existence of the BO phase 
in our simulations. One observation is the "phase sep- 
aration" or "flip-flop' between left and right BO phases 
as a function of imaginary time in the QMC simulations. 
The other is the staggered behavior of the BO correlation 
functions for which the numerical data out to large dis- 
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tance supports power law behavior in the non-ionic case 
and long range order in the ionic case. 

Other work on related systems also has identified BO 
phases. Recent work by Nakamural^ in the extended 
Hubbard model has found a rich phase diagram in which 
there are 2 transitions from a BI BO and BO MI 
regime. The extended Hubbard model differs from the 
ionic model studied here in that there is an additional 
next nearest neighbor coulomb potential V and no ionic 
potential A. Nakamura identifies the first transition as 
belonging to the gaussian universality class and the later 
as a Kosterlitz-Thouless transition. The BO phase ob- 
served by Nakamura exists at all V down to the usual 
Hubbard model {V = 0); where both the gaussian and 
the KT transitions coexist. We can imagine the Uc at 
which the BI BO transition takes place increasing 
concurrently as the ionic potential is increased from zero. 

Let us now consider why the BO phase was not found 
in previous studies that used exact diagonalization Lanc- 
zos techniques to treat small finite systemsElS. There 
are two reasons why these studies did not find the BO 
phase. In the BO phase the energy can parameterized by 
the bond order parameter (Eq |^) that develops a bimodal 
distribution with minima at ±_B. The true ground state 
is a linear combination of these 2 degenerate BO states, 

*o = cos(6')«'+ + sin(6i)1'_, (29) 

which of course has no net bond order. The situation 
is similar in many aspects to a ferromagnet; it is only in 
the thermodynamic limit that one or the other of the two 
states is the true ground state with long range order. For 
finite systems existence of the BO state can be inferred 
from correlation functions; however, to our knowledge 
this has not been done in other work. 

A second reason that the BO states have not been ob- 
served may be that there is no bimodal distribution for 
the small cells studied by exact diagonalization. We have 
addressed this issue using QMC by measuring the aver- 
age bond order on distorted lattices of 14 < L < 62 sites 



and extrapolating to the centrosymmetric limit. At this 
point in the phase diagram (A = 0.0716, U = 2.4) lat- 
tices with less than 50 sites do not exhibit the BO phase 
and only upon working with larger supercells does QMC 
detect the phase separation of the two BO states. Conse- 
quently, exact diagonalization methods are not currently 
feasible in such cases since they scale exponentially with 
L. 

Recent work using the density matrix renormalization 
group (DMRG) method has reported results for charge 
Ac and spin As gaps in these models. This approach 
should enable one to distinguish the phases since (i) Ac = 
As ^ in the BI phase, (ii) Ac ^ ^ in the BO 
phase, and (iii) Ac ^ but A^ = in the MI phase. It 
was found to be very difficult and to require extremely 
large cells to determine spin gaps in the BO/MI phases, 
and the two reports came to opposite conclusions on the 
existence of the BO phase. In our QMC calculations 
we have also determined the charge and spin gaps. Our 
estimates of the charge gap are in qualitative agreement 
with other works; however, the spin gap is very small in 
all cases except in the BI regime and statistical noise does 
not permit accurate determination of such small gaps in 
QMC. 

Both DMRG calculations find the spin gap to vanish, 
i.e., the MI phase to be the ground state for large U. We 
have no direct explanation of this difference: it may be 
that our procedure is not sufficiently accurate to deter- 
mine the BO-MI transition, which is the most difficult 
part of the present work. On the other hand, it may 
be that the DMRG calculations on finite cells with open 
boundary conditions may have difficulties: the surface 
effects break the symmetry of the problem which may 
lead to extremely problematic size effects and potential 
errors. In any case, we are very confident that our work 
establishes that the BO state is either the ground state 
or very close to the ground state in energy; this is clear 
from our tests on the ordinary Hubbard model shown in 
Fig. 1. 
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VIII. CONCLUSIONS 

We have studied the phase diagram of an ideahzed 
dielectric, the ID ionic Hubbard model proposed by 
Nagosail and Egamill. This model undergoes a phase 
transition as a function of the on-site interaction U, 
which has been a source of controversy. The only previ- 
ous quantitative studiesEBil concluded that at a critical 
U there is an abrupt " topological" transition from a band 
insulator to a Mott insulator with no broken symmetry 
or long range order in either phase. The signature of the 
transition was found to be an abrupt change of 1/2 in the 
polarization at which the effective charge diverged signi- 
fying the delocalization of the electron stateJllS. Re- 
cently, however, there has been a predictionll that this 
model would exhibit two quantum phase transitions: the 
first signifying a change of state from a band insulator 
to a broken symmetry phase with long range alternat- 
ing bond order, and the second a transition to the Mott 
insulator. 

We have studied this model using quantum Monte 
Carlo methods which allow the simulation of much larger 
systems than studied by exact diagonalizationEi@. To 
our knowledge, this is the first application of QMC to 
determine the polarization and localization of an elec- 
tronic system. We evaluate the expectation values of the 
bond-order, polarization and localization using the ex- 
pressions Eqs H, 1^ and |^. It is found that upon crossing 
a critical value Uc a change of phase occurs from a band 
insulating to bond-ordered state. The bond order de- 
velops continuously (See Fig ||) as a function oi U — Uc 
and since the inversion symmetry is broken, the polar- 
ization also varies continuously, unlike the results of the 
small cell exact-diagonalization calculations .111 The crit- 
ical behavior is uniquely determined by fitting the bond 
order and polarization to a scaling function near the crit- 
ical regime. We find an exponent near 1/2, which differs 
from that for the Ising class proposed in Refi; however, 
it may be that we are outside of the regime in which the 



scaling belongs to the appropriate universality class. In 
addition, we found that there is a metallic point at Uc 
where the system is metallic. At this point the charge gap 
must vanish which we have found in pure ground state 
calculations by determining the fluctuations of the po- 
larization. The calculations determine quantitatively the 
localization lengthllUli, which diverges at the transition. 

An important part of the present QMC work is that we 
use a forward projection scheme which allows exact esti- 
mates, in principle, of any operator, including ones such 
as the polarization (or center of mass position operator) 
that do not commute with the Hamiltonian. Further- 
more the nodes of this ID model are known exactly, so 
the QMC method is in principle exact. In order to con- 
firm the existence of the bond ordered state, we carried 
out calculations on dimerized lattices (to ± 6) whose in- 
version symmetry is explicitly broken, and let S become 
small. QMC allows us to work with large enough lat- 
tices so as to study systems with levels of dimerization 
much smaller than previously feasible0~00. We find 
good agreement with previous results obtained from the 
Heisenberg spin model that predict electronic correlation 
enhances the instability to bond ordering. In addi- 
tion, we can see from the simulations of the symmetric 
lattice that the system is alternating between the two 
degenerate states of bond-order (see Fig^). 

We have searched for the proposed transition to a Mott 
insulating state, but we have not observed such a tran- 
sition from the bond ordered regime even for very large 
U or very small A. Even the smallest value of A consid- 
ered in this study (A/to = 1/14 < U/to = 2.4) is suf- 
ficient to cause the ionic Hubbard model to be unstable 
to bond ordering, although there is no broken symmetry 
in the usual Hubbard model (A = 0), neither in the ex- 
act solutionii nor in our results. Thus our results show 
that the instability to dimerization is even stronger in 
the ionic model than that known previously for the ordi- 
nary non-ionic Hubbard modeLElFilEill^i. Furthermore, 
for the centrosymmetric lattice {S = 0), calculations of 
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correlation functions and observations of "flip-flop" be- 
tween left and right bond-ordered states in the QMC sim- 
ulations provides further evidence for the stability of the 
bond-ordered state. 

Among the interesting consequences of the stability of 
the BO state is the existence of fractional charges. BE 
For the case of a dimerized or bond-ordered state, the 
charge is an irrational fraction the value of which depends 
upon the value of aEIEI. 

Finally, these results imply that if dimerization is al- 
lowed (which is always the case in real materials since 
the atoms can always dimerize if it lowers the energy) 
then the symmetric Mott state is never stable and the 
only phase transition is from the symmetric BI to the 
dimerized BO state. This is experimentally confirmed 
by Torrance et al@; where upon increasing the electronic 
interaction a BI BO transition takes place. 



TABLE I. Fitting parameters for polarization and bond 
order. The quantities in parenthesis are the error in the last 
decimal place. 
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